Shape and dynamics of seepage erosion in a horizontal granular bed 
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We investigate erosion patterns observed in a horizontal granular bed resulting from seepage of 
water motivated by observation of beach rills and channel growth in larger scale landforms. Our 
experimental apparatus consists of a wide rectangular box filled with glass beads with a narrow 
opening in one of the side walls from which eroded grains can exit. Quantitative data on the shape 
of the pattern and erosion dynamics are obtained with a laser-aided topography technique. We 
show that the spatial distribution of the source of groundwater can significantly impact the shape of 
observed patterns. An elongated channel is observed to grow upstream when groundwater is injected 
at a boundary adjacent to a reservoir held at constant height. An amphitheater (semi-circular) 
shape is observed when uniform rainfall infiltrates the granular bed to maintain a water table. 
Bifurcations are observed as the channels grow in response to the ground water. We further find 
that the channels grow by discrete avalanches as the height of the granular bed is increased above the 
capillary rise, causing the deeper channels to have rougher fronts. The spatio-temporal distribution 
of avalanches increase with bed height when partial saturation of the bed leads to cohesion between 
grains. However, the overall shape of the channels is observed to remain unaffected indicating that 
seepage erosion is robust to perturbation of the erosion front. 

PACS numbers: 45.70.Qj, 47.56. +r, 47.57.Gc, 92.40.Pb 



I. INTRODUCTION 



The shape of rivers derives from their complex inter- 
actions with landscapes. The flow of water erodes the 
land which in turn modifies the flow producing meanders 
and networks of rivers [TH3] . Of significant interest is the 
case of seepage erosion or sapping [4H6] in which a flow 
of water occurs inside a permeable and erodible layer 
of sediment above an impermeable layer of a different 
composition. Once the water inside this porous medium 
emerges at the free surface, producing a spring, the flow 
removes grains from the surface by erosion, progressively 
digging a deeper channel, which in turn can draw more 
water, inducing the growth of a river. Seepage erosion is 
said to shape many examples of valleys, canyons and river 
networks and assumed to produce amphitheater-headed 
valleys HG3IZHS1, although recent field reports [101 ITT] 
suggest that some of these examples can be attributed to 
overland flow as well. At smaller and more rapid scales, 
seepage erosion has been considered in the formation of 
channels on the beach during outgoing tide [5j IT2UT4] . 
where capillary cohesion between grains can be also im- 
portant. The spatial distribution of the groundwater flow 
when brought from a far away reservoir as opposed to be- 
ing fed by local precipitation can be crucial to the shape 
of the channels. 
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While the evolution of a river is highly non-linear, the 
geological variability of the land, initial topography, and 
changing climate can also make it difficult to identify 
and extract common features important to their initial 
formation and growth. Laboratory experiments can help 
unravel fundamental physical phenomena involved in ero- 
sion and the development of channels. Indeed, seep- 
age erosion, channel formation and bifurcation has been 
demonstrated with experiments by Howard and collab- 
orators using an inclined granular bed with subsurface 
water flow injected at a boundary [I5j [16], and by a lo- 
calized and continuous precipitation (mist) [17]. Further 
work by our group has shown the effect of slope of the ini- 
tial bed and the height of the water table on the spacing 
of channels observed [14]. We have also used the Shields 
number, which captures the ratio of the viscous drag and 
gravitational forces, to explain the onset of channeliza- 
tion and slumping of the bed, and its dependence on 
driving conditions [18]. Quantitative topography mea- 
surements of the erosion patterns were also used to pro- 
pose a dynamical model for channelization [T9] . 

Recent experiments by Izumi and collaborators [2Q| [2Tj 
have also shown channel bifurcation by using coarse sand 
and a weak slope [20 . The authors attribute this occur- 
rence to the important role of the geometry of ground- 
water flow and the "resistibility of sediment material 
to slope failure." Using a similar experimental appara- 
tus, but with overland flow on cohesive soil [22 , the 
same group also observed channel bifurcations. Finally, 
complex channels produced by a combination of seepage 
and surface runoff flows have been also investigated [23] . 
Thus, it remains unclear if overall channel shape can be 
used to extract if seepage or overland flow is dominant 
in a given situation and a detailed investigation of both 
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situations is necessary to understand similarities and dif- 
ferences. 

In this paper, we focus on a situation where the gran- 
ular bed is initially flat and the water table is fed by 
uniform rain to further connect to examples of seepage 
channels given in the field, which appear to grow in flat 
landscapes [HE]. To our knowledge, seepage erosion of a 
flat layer of sand has been never studied experimentally. 
Although this case may appear simpler than the inclined 
case, it is in fact more complex because the bed slope be- 
comes a free parameter evolving during the growth of the 
erosion pattern and influencing system dynamics. Fur- 
ther, we study here only the case of erosion driven by 
seepage, and eliminate overland flow to simplify the prob- 
lem. In particular, we seek to understand the influence 
of the rainfall and the resulting groundwater distribution 
on the shape and dynamics of resulting erosion patterns. 

The paper is organized as follows. We first introduce 
the essential feature of seepage erosion in Sec. [TTJ and 
then we present the experimental apparatus used in our 
investigation in Sec. |III| We then compare in Sec. IV the 
obtained erosion patterns when seepage flow is driven by 
rainfall and by a flow inject at the boundary to simulate 
groundwater dominated by precipitation outside the box 
containing the granular bed in the upstream direction. 
Then, the effect of depth of the eroded bed is also in- 
vestigated, enabling us to demonstrate that the erosion 
driving mechanisms are robust in presence of surface per- 
turbations. From these measurements, we discuss the in- 
fluence of rainfall on the shapes of patterns and growth 
dynamics produced by seepage erosion in Sec. V. 



II. PHYSICS OF SEEPAGE EROSION 

A schematic of a model seepage erosion geometry 
is shown on Fig. [I] We consider a horizontal "step" 
of porous erodible medium, placed on an impermeable 
and non-erodible layer. Water moves inside the porous 
medium fed by uniform rain above with a rate which is 
lower than the infiltration rate to prevent overland flow. 
Further, the water can be also fed through the boundary 
from the upstream direction due to precipitation further 
upstream or a reservoir. The water table corresponds 
to the height inside the bed where water saturates the 
medium and is tilted depending on the flow present, the 
porosity of the bed, and the boundary conditions. A 
spring appears where the water table intersect the bed 
surface, leading to surface flow. The shape of the bank 
above the surface flow is set by its mechanical stability 
and conservation of mass. In case of unconsolidated bed 
composed of sedimentary sand, the slope of the bank is 
usually given by its angle of repose [19] . The surface flow 
of water erodes the granular bed by removing grains at 
the bottom of the stream depending on flow rate [24] [25] , 
which in turn destabilizes the bank and leads to growth 
of a channel in the upstream direction. While these over- 
all features are well accepted, a detailed understanding of 
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FIG. 1: (Color online) A schematic cross-sectional view of the 
geometry of seepage erosion. Water moves inside the bed from 
upstream (right) to downstream (left) direction. The water 
table is the surface separating medium saturated by water 
from the unsaturated medium, and intersects the granular 
surface where the water flow emerges at the surface. Grains 
can be removed by water flowing near the surface leading to 
headward growth in the upstream direction. 



the angle at which flow emerges at the surface, its effect 
on the erosion rate, and development of the erosion front 
in response to surface perturbations is far from clear. 

Away from the surface, the flow inside the porous 
medium is modeled by Darcy's law with appropriate 
boundary conditions. Further, when the vertical com- 
ponent of the flux is small relative to horizontal com- 
ponents, the flow in the porous medium can be approx- 
imated as two dimensional, which is called the Dupuit 
approximation [6 . Therefore the three dimensional prob- 
lem can be reduced to two dimensions and the seep- 
age flow can be obtained from the water table elevation 
h(x,y). In the absence of rain, it can be demonstrated 
that h 2 follows the Laplace equation [6] 

V 2 h 2 = 0. (1) 

with appropriate boundary conditions. In the presence 
of homogeneous rain, a source term transforms this equa- 
tion into a Poisson equation 

V 2 h 2 = ~, (2) 

where, P is the rainfall rate per unit surface area pro- 
jected on the horizontal plane, and k the hydraulic con- 
ductivity of the porous medium. Because of the impor- 
tant differences in the mathematical properties of these 
two equations, the shape of the resulting groundwater 
flow and thus the resulting erosion dynamics are expected 
to depend strongly on the relative importance of local 
rainfall compared to upstream flow. 

In order to illustrate the differences in the groundwa- 
ter flow near a channel under the two different driving 
conditions, we obtain the shape of the water table by 
simulating Eq.(l) and Eq.(2) with appropriate boundary 
conditions [see Fig.[2](a,b)]. We then find the groundwa- 
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FIG. 2: (Color online) (a) Simulation of Eq. [I] in a domain 
in the presence of a channel. Water is injected on the top 
boundary (h = 1) and flows through the bottom boundary 
with a channel (h = 0). Contours correspond to h 2 and the 
arrows represent the streamlines of the groundwater flow, (b) 
Simulation of Eq. [2] in presence of a channel. Water is rain- 
ing homogeneously on the domain and flows through the top 
and the bottom boundaries, (c) Groundwater flux along the 
arclength of the channel length corresponding to the simu- 
lation for the Laplace equation. Maximal flux occurs at the 
tip of the channel, (d) Groundwater flux along the channel 
arclength corresponding to the Poisson equation shows two 
symmetrically placed maxima. 




FIG. 3: (Color online) Schematic of the experimental appa- 
ratus. A Plexiglas box which is open at the top is filled with 
glass beads with diameter 0.5 mm. The front and back side 
of the granular bed are in contact with water reservoirs main- 
tained at heights Hi and H2 above a boundary which allow 
water to pass but not the grains. An opening into the front 
reservoir allows the eroded grains to exit the box. The bed 
can be exposed to uniform rain produced by an array of three 
nozzles (not represented). Channel depth corresponds to the 
vertical elevation difference from the initial flat bed surface 
height to the bed surface. Topography is measured inside the 
observation window using a laser scanning technique. Dimen- 
sions are given in the text. 



ter flow rate per unit width using Darcy's law: 

Q(x,y)Tt = -^h 2 (x,y), (3) 

where, it is the unit vector parallel to ^h 2 (x,y). We 
note that in the Laplace case, flow entering the channel is 
mostly localized near its tip [see Fig. [2^c)]. Whereas, the 
groundwater enters the channel almost uniformly from 
all directions in the Poisson case, and have two symmet- 
rically placed maxima [see Fig. [2^d)] . Further, we have 
shown in the past [26] that the local erosion velocity at 
a channel front increases linearly with the groundwater 
flow entering at that point. Therefore, when the ground- 
water flow is driven by an upstream flow (Laplace case), 
we expect that the channel grows preferentially at its 
tip. By contrast, water is collected more uniformly on 
the channel boundary in presence of rain and the max- 
ima of groundwater flux are located on the sides of the 
channel. Thus, higher erosion rates can be expected near 
these regions leading to channel splitting. 

Although the thickness of the permeable bed does not 
affect the seepage flow itself, it does impact the mechan- 
ical stability of the bank and the overall channel shape 
that can develop. In the inclined case, the initial bed 
slope is the parameter determining the mechanical sta- 
bility of the sediment [14 . For the flat bed considered 
here, an analogous parameter may correspond to the ra- 
tio between the largest depth of granular material above 
the water table and the maximum length of the channel 
and corresponds to the minimum slope inside the chan- 



nel. Moreover, for small-scale erosion patterns like beach 
rills, mechanical stability should be strongly affected by 
capillary cohesion [27]. Above the water table the un- 
saturated granular medium can become indeed cohesive 
due to capillary bonds between particles [28, 29 , which 
are present if the sand is wet. By contrast, the satu- 
rated granular medium below the water table remains 
cohesionless. This vertical inhomogeneity of the bed can 
influence how the material from the bank erodes into the 
channel [30], before it is removed by the water flow. The 
occurrence and distribution of these bank collapses de- 
termines the overall shape of channel. 



III. EXPERIMENTAL SETUP 

The experimental apparatus consists of a Plexiglas box 
(120 cm wide and 120 cm long) which is filled with cohe- 
sionless glass beads with diameter 0.5 mm as shown in 
Fig. [3] The hydraulic conductivity k for this porous ma- 
terial was measured to be 3mms _1 [13] , and the height 
of capillary rise of water d cap is measured to be 3.2 cm. 
Practically, d cap gives a minimum depth of the outlet 
doutiet to avoid surface flow in the experiments because 
it gives the height of the water table above where the flow 
emerges from the granular bed. We label the direction 
upstream from the outlet as the y axis, its perpendicular 
in the horizontal as the x axis, and the vertical height 
as z. A metal blade on a track is used to set the initial 
depth of the bed at d san d = 17.5 cm. The bed is also 
confined (see Fig. |3| in the y direction between two walls 
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FIG. 4: Images of the four kinds of channels observed in the experiments, (a) An elongated channel fed by a flow from a 
reservoir at the back of the bed after t = 916 min with d ou tiet = 4 cm, Sh = 5 cm and no rainfall (P = 0). (b) A wider 
bifurcated channel fed by rain (Pi ~ 1.04 x 10 -4 ms -1 ) is observed with d ou tiet — 4cm (Sh = 5cm, t = 595 min). (c) An 
elongated deeper channel with steeper banks is observed with d ou tiet = 8.0cm, (P = 0, Sh = 6.8cm, t = 311 min). (d) The 
channel is deeper with doutiet — 8.0cm and wider with rain (P = Pi, Sh = 2cm, t = 240 min). 



separated by distance 66.5 cm, which have a mesh that 
allow water to pass but not the grains. The boundary in 
the front with the outlet for sand is designed to be per- 
meable to water in order to create a subsurface flow from 
the back to the front which is unperturbed by the outlet 
at large distances along the positive and negative x di- 
rections. The outlet consists of a 6 cm wide rectangular 
slot as shown in Fig. [3j with a depth d out i e t which can 
be varied continuously between 11 cm and 3.5 cm. Water 
reservoirs with heights Hi and H<i on the other side of 
the front and back permeable walls can be adjusted using 
a water pump to induce a subsurface flow from the back 
to the front of the bed. 

We also use commercial nozzles which produce small 
water droplets in a full cone with a 90 degree aperture to 
mimic rainfall on the bed itself. To produce a homoge- 
neous and controlled rainfall rate poses some experimen- 
tal challenges. After experimentation, we choose to use 
three nozzles separated by 15 cm and placed above the 
surface at a distance H ra i n = 46 cm off the sand surface. 
The y coordinate of the central nozzle is L rain = 35 cm. 
We obtained a fairly homogeneous rain in the central re- 
gion of the bed and within 30% overall decrease towards 
the edges. The rain itself can produce a supplemen- 
tary erosive process of the granular bed. For example, 
rainfall on a cohesive material can create complex land- 
scapes [31]. In our experiments, the impact of the water 
drops act to only smoothen spatial variations of the gran- 
ular surface and does not modify the seepage flow [32] . 
Experiments with rain are conducted with a total flow 
rate of 69 cm 3 s _1 with 10% fluctuations. The corre- 
sponding rainfall rate Pi ~ 1.04 x 10 -4 ms _1 is obtained 



as the ratio of rain per unit time and the surface area in 
a central zone. To avoid any excess water at the bound- 
aries of the box, a gutter is also placed right next to it to 
remove water which hits the side boundary. Finally due 
to the presence of this gutter and the blade system, the 
actual length in the y direction over which experiments 
are conducted is L = 55 cm. 



A number of parameters can be varied in our experi- 
ments. These include the size of the outlet d out i e t, which 
controls the eroded depth by setting its maximal pos- 
sible value, the amplitude of the upstream flow by set- 
ting the difference between H^ and Hi , Sh = Hi — Hi . 
In order to simplify the experimental situation, we set 
Hi to the same level as the bottom of the outlet i.e. 
Hi = d san d — d ou tiet' Further, presence or absence of 
rainfall is simulated by turning on and off water circula- 
tion to the nozzles. The amplitude of the upstream flow 
rate due to Sh is proportional in the framework of Dupuit 

2 2 

approximation with the parameter U = ^ [6], 



L 2 

where L is the length of the bed. This corresponds to a 
dimensionless seepage velocity due to the upstream flow 
alone. We note that for the Dupuit approximation to be 
valid, we need U <C 1. Finally the dimensionless rainfall 
2P 

rate is — — and corresponds to or 0.069 depending on 
k 

the absence and presence of local rain in our experiment. 



We try to estimate the flow rate at the position of the 
outlet (x = 0, y = 0), using the Darcy law, Eq. [3] . With 
the approximation that the seepage flow depends mainly 
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on y, we obtain: 

with a difference in height of the front and back reservoir 
and without uniform rain. With only uniform rain, we 
obtain: 

Qr == PL- 

where, L is the length of the bed in the y direction, and W 
is the total width of the sand box. For a significant height 
difference Hi = 9.5 cm and H2 = 16.3 cm, the flow rate 
can be estimated to be Q w ~5 x 10 -6 m 2 s _1 . Whereas 
the flow rate due to rain without height difference is Q r ~ 
1 x 10 _4 m 2 s _1 . Therefore, the seepage flow rate is at 
least one order of magnitude larger due to rain and thus 
dominant in our experiment. 

The surface is observed with a CMOS camera located 
approximately 2 m above the granular bed. For quan- 
titative measurements of channel shape, we use a laser 
scanning technique [19] to reconstruct the topography 
inside an observation window (-40 cm < x < 40 cm and 
5 cm < y < 40 cm) indicated in Fig. [3] A laser sheet is 
swept over the bed surface along the y axis and the illu- 
minated crossection of the bed at that location is imaged 
and used to reconstruct the topographic map of the bed 
after appropriate calibration. The time when erosion is 
observed to start is chosen as time t = 0, and is studied 
over several hours with successive scans. Experimental 
analysis is terminated when channel size exceeds the size 
of observation window, which also avoids direct effect of 
the boundaries on the growth of the channels. We also 
confirmed that the presence of water flow at the surface 
few millimeter deep does not affect the topography mea- 
surements. 



IV. OBSERVATIONS 

For a sufficiently high seepage rate, erosion is observed 
near the outlet with a sediment flux rate which is related 
to the seepage flow rate. The erosion front can stop after 
progressing in the upstream direction or continue till the 
boundary is reached depending on the applied conditions. 
We call the resulting patterns channels because they are 
formed as a result of erosion and a stream is present at 
the bottom of the eroded bed, even though the ratio of 
the length and the width may be comparable and not as 
high as in case of channels in the conventional sense. The 
grains removed through the outlet do not play a further 
role in the erosion back stream. For sufficiently large 5h, 
or in cases where d out i e t < d capi the entire granular bed 
becomes saturated and a surface flow occurs over the 
bed which is different from seepage erosion [TJ |2], and 
therefore not discussed further here. 

Because of the bed geometry and the difficulty in ac- 
tually measuring the water flow on and near the surface, 



it is not possible to compute the relevant Shields num- 
ber. We notice nevertheless that in our experiments, sur- 
face flows cover a large area of channels with a depth of 
a few millimeters which is at least one order of magni- 
tude smaller than the channel depth. Consequently, we 
focus our study on shape evolution and growth of chan- 
nels as a function of the type of seepage flow and as a 
function of channel depth controlled by the height of the 
outlet d ou tiet- A variety of channel shapes and growth 
velocities can be obtained using the reservoir or uniform 
rainfall and various outlet depths. Figure [4] shows im- 
ages of four different examples obtained by using various 
combination of outlet depth and sources of seepage flow. 



A. Effect of the spatial distribution of seepage 
water source on channelization 

We begin by discussing the examples for small val- 
ues of d ou tiet where the Dupuit approximation is valid, 
and where the channel depth is just above the threshold 
needed to prevent overland flow. The map of a resulting 
channel depth d(x,y) relative to the initial surface ob- 
tained with laser aided topography after time t = 677 min 
is shown in Fig. |5|a) for d ou tiet — 4.0 cm, Sh = 5 cm, 
and P = 0. A smooth finger shaped channel is ob- 
tained with more or less uniform internal slope. To il- 
lustrate the evolution of the shape, we plot the contour 
d{x, y) = 1 cm corresponding to the external boundary 
of the channel at various times in Fig. |5|c). The growth 
appears directed to the back of the bed and the channel 
width remains constant. By contrast a significantly wider 
channel is obtained when groundwater is fed by uniform 
rain (see Fig. [5J3). The evolution of the corresponding 
d(x, y) = 1 cm contour is shown in Fig. [5^d). We observe 
that the growth is isotropic initially, but appears to then 
grow towards the top left and right corners as the chan- 
nel starts to grow after a distance of about 15 cm from 
the outlet. 

As illustrated in Fig. [2|b,d), we expect that the 
groundwater flux entering into the channel has two max- 
ima on each channel sides in the presence of homogeneous 
rain. As erosion rate should be higher at these points, a 
channel splitting could be so induced. This bifurcation 
could be the first step toward creation of a river net- 
work, such as the seepage networks of the Florida Pan- 
handle [8j|9]. However, the typical channel width (which 
is set by the flow rate and the size of grains) [33] is large 
compared to the system size to observe further channel 
splitting. By repeating experiments under similar condi- 
tions, we find that the bifurcation always occurs but in 
general not symmetrically and can sometimes take more 
complex shapes. This variability may indicate well an 
instability in addition to a bifurcation driven by ground- 
water flow. Finally, the small scale topography inside 
the channel visible on the pictures Fig. [4] (a-b), appears 
more complex in the case where P = 0. The surface flow 
is greater when P = Pi , further eroding the channel bot- 
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FIG. 5: (Color online) Comparison of channels in a shallow channel regime where the capillary rise is of order of the channel 
depth, (a) Depth map of the channel observed after t = 677 min with Sh = 5 cm, P = and d ou tiet = 4 cm. Depth is indicated 
with a scale in centimeters measured from the initial flat bed. (b) Depth map of channel observed after t = 300 min with 
Pi - 1.04 x 10" 4 ms"\ S H = 0.9 cm and d ou tiet = 4 cm. (c, d) The corresponding channel shape evolution for the depth 
d(x, y) = 1 cm contour increase in length with time. 
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FIG. 6: (Color online) (a,b) Channel length and eroded vol- 
ume evolution in absence of rain P = with Sh = 5 cm. (c,d) 
Channel length and eroded volume evolution in presence of 
rain P — Pi with Sh — 0.9 cm. 



torn, making it nearly flat in our experiments. Falling 
droplets also appear to smooth small surface perturba- 
tions [32 and the channel bank appears more regular 
with rain. 

To characterize the channel dynamics, we use the chan- 
nel length L c h defined as the distance from the outlet to a 
point along the y-ax.is where the depth is above 1 cm i.e. 
d(0,L c h) = 1cm. Fig. [6^a,c) shows the channel length 
as a function of time for the case with imposed upstream 
flow and rainfall, respective. One can observe that in 
case of upstream flow, the graph has an upward curva- 
ture indicating that the channel velocity increases as it 
grows. This can be understood from the fact that the 
channel tip approaches the boundary where water is in- 



jected. Because a significant amount of water flows be- 
low the channel surface, the closer the channel is to the 
source of upstream flow, the greater is the relative im- 
portance of surface flow which results in greater erosion. 
On the other hand, the curve for the rainfall case has a 
downward curvature. This occurs because a significant 
fraction of the rain starts to fall inside the channel de- 
creasing the seepage contribution which in turn reduced 
the growth of the channel. Because the choice of the di- 
rection is somewhat arbitrary and may be influenced by 
the splitting of the channel, we also calculate a measure 
of the scale of the channel using the volume of the grains 
eroded in the channel as a function of time. The data 
were obtained by integrating the measured depth field 
and are plotted for the two cases in Fig. gM). The 
trends shown by the channel length graphs are reflected 
in these graphs as well. One sees that the velocity of the 
length grows less slowly compared with the erosion rate, 
indicating that the channel with rain grows wider with 
time, indicating a bifurcation. 



B. Effect of cohesion on channelization 

We next investigate how channelization occurs when 
the channel depth is large compared to capillary rise, 
and a cohesive layer is present above the water table be- 
cause of capillary bonds between grains. Due to seepage 
flow below the unsaturated bed which undermines the 
bed, the bank becomes mechanically unstable and col- 
lapses [30] . We first describe the global effect of increas- 
ing the depth on seepage erosion before characterizing 
the nature of the bank avalanches. 
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FIG. 7: (Color online) Comparison of channel shape in the regime where d ou tiet > d cap and significant layer of unsaturated 
cohesive bed layer is present, (a) Depth map of channel observed after t = 306 min with d ou tiet = 8.0 cm, Sh = 6.8 cm and 
P = 0. (b) Depth map of channel observed after t — 225 min with d ou tiet = 8.0 cm, Sh = 2.0 cm and P — P\. (c, d) The 
corresponding channel shape evolution for the depth d(x,y) = 1cm contour. The contours increase in length with time. 
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FIG. 8: (Color online) Comparison of channelization dynam- 
ics where d ou tiet > d cap when the seepage erosion is gener- 
ated by an upstream flow and by rain. (a,b) Channel length 
and eroded volume evolution for P = with d ou tiet — 8 cm, 
Sh = 6.8 cm. (c,d) Channel length and eroded volume evolu- 
tion for P — P\ with d ou tiet — 8 cm, Sh — 2.0 cm. 



Figure [T^a) shows an example of a channel formed with 
seepage flow derived only from the reservoir (d out i e t — 
8.0 cm and Sh = 6.8 cm). The shape of the channel re- 
mains narrow and it grows towards the back of the bed 
with a constant width (in Fig. [7jc)). Thus, the overall 
shape is similar to that observed with a shallower out- 
let. However, there are differences in how the erosion 
front advances. For example, if one compares the con- 
tours for t = 157 min and t = 255 min, the two contours 
are identical except around 10 cm. This local advance of 
the contour can be observed at other times as well. By 
contrast, Fig. [7^b,d) shows examples of channels formed 
with uniform rainfall on the bed with d out i e t — 8.0 cm, 



and Sh — 2.0 cm. In this case, the channels appear more 
or less circular like amphitheater-headed valleys observed 
in nature in agreement with the notion of an isotropic 
growth in presence of homogeneous rain. However, as 
the erosion front becomes comparable to the size of the 
bed, the growth appears more towards the top left and 
right corners as in the case for shallow channels shown 
in Fig (5^c). In contrast with shallow channel discussed 
earlier, a bifurcation is not observed possibly because the 
deeper channels draw more water from the reservoir and 
still continue to grow towards it. 

Next, plotting the channel length for these two cases 
(Fig. [8^a,c)), we find that the channel length grows in 
discrete steps with a decreasing amplitude and time in- 
terval between successive avalanches when seepage flow 
is imposed with a reservoir, whereas the avalanches have 
a larger amplitude and are less frequent in cases with 
uniform rain. The temporal evolution of volume (Fig. [8] 
(b,d)) shows similar trends, but without the steps. Al- 
though material breaks from the bank of the channel, it 
falls to the bottom of channel and is then removed by the 
surface flow. Consequently, although rapid avalanches 
change the local channel shape, materials leave the sys- 
tem at a lower steady pace due to surface flow. Fur- 
ther, if we compare volume evolution in the two cases 
shown where flow is imposed from a reservoir (Fig. |6|b) 
and Fig. [8^b)), we notice that the curves have similar 
shapes. Comparing the cases with rain (Fig. |6^d) and 
Fig. |8jd)), we find that the volume of removed sediment 
is roughly linear for deeper channels, and therefore the 
erosion rate is constant during the experiment. Thus, 
slowing of erosion velocity is not observed and it appears 
that the reduction of the surface exposed to rain does not 
influence erosion rates for deeper channels. This observa- 
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FIG. 9: (Color online) (a) Sequence of images illustrate a 
channel bank collapse. A crack initially appears and breaks a 
block, which falls inside the channel over the time scale of a 
minute, (b) The spatio-temporal development of avalanches 
as a function of the angle and time (b) in case of flow from 
reservoir and (c) in case with local rain (the experimental 
parameters are the same as in Fig. The colorscale is in 
arbitrary units and corresponds to the temporal derivative of 
the distance from the contour at d(x, y) — 1 cm and origin. 



Because channel length shown in Fig [8^a,c) detects 
avalanches only on the central axis, we now use mea- 
surement of the full contour defined at an eroded depth 
d(x, y) — 1 cm. For each point of this contour M, we 
compute the distance from origin (1m as a function of the 
angle to the ?/-axis. The time derivative of (1m as a func- 
tion of position and time is plotted in Fig.[9jb,c). In case 
of flow from reservoir, avalanches are homogeneously dis- 
tributed at the beginning and then become increasingly 
localized in the central direction with decreasing ampli- 
tude and increasing frequency. In contrast, avalanches 
are more regular and of uniform amplitude under homo- 
geneous rain. Collapse events start at the center and 
then propagate to the edges over a timescale of about 
ten minutes. 

Although the dynamics are very different in the two 
cases, we note that the time between two avalanches de- 
creases as the channel grows when the flow comes from 
a reservoir. In both cases, it appears that the time in- 
creases when the amplitude of previous avalanches in- 
creases. This is because the channel has to advance fur- 
ther to destabilize the bank after a large avalanche similar 
to observations in horizontally rotated cylinders [33] . In 
the reservoir derived flow, the internal slope inside the 
channel increases progressively with distance to the ori- 
gin. Therefore, the vertical size of a block which can fall 
decreases during channel growth and avalanches become 
smaller and more frequent. In spite of the intermittent 
nature of local growth, the overall shape is determined 
by the distribution of the source of seepage flow over the 
duration of our experiments. 



V. ANALYSIS OF CHANNEL SHAPES 



tion further suggests that upstream flow gives significant 
contribution to seepage erosion in case of deeper chan- 
nels as they approach the reservoir, and suppresses the 
bifurcation of the channel. 

An example of a bank collapse in this cohesive regime is 
shown in Fig. |9|a) and appears similar to those reported 
previously [18 , 30 . A block of grains of typical size 10 cm 
- held together by capillary bonds - is observed to break 
off with a vertical crack and falls inside the channels. Suc- 
cessive such events produces an irregularly shaped bank 
with a steep wall as can be observed in Fig. [7] The break- 
ing of a block of material can be understood as follows. 
As the flow in the saturated region of the bed erodes 
grains, the cohesive strength of the unsaturated region 
on top (see Fig. TJ leads the surface to be stable till a 
threshold is reached. This threshold leads the bank to 
fail in blocks whose size is given by a balance of the vis- 
cous drag below and gravitation which tend to destabilize 
the bank, and the capillary forces which tend to stabilize 
it above the angle of repose of the dry grains [27] . 



In order to examine the effect of the experimental con- 
trol parameters on the erosion rate and the aspect ratio of 



the channel, we plot these quantities in Fig. 10 We find 
that U has to exceed a threshold « 0.04 for erosion when 
P = 0, and erosion occurs for even U — when P = P\. 
The erosion rate r s plotted in Fig. [lO^a) is computed as 
the time derivative of the channel volume obtained from 
the topographic maps. In order to compare different ex- 
periments, we show r s corresponding to channel lengths 
equal to 20 cm. For experiments without rain, the av- 
erage erosion rate < r s >= 3.1cm 3 min _1 and the cor- 
responding standard deviation is o Ts = 1.9cm 3 min _1 . 
With rain these values become < r s >= 30cm 3 min _1 
and a Ts = 20cm 3 min _1 . Therefore the erosion rate is 
one order of magnitude higher in the presence of rain for 
the given rainfall rate. But the high value of the stan- 
dard deviation shows that other parameters play a role 
as well. The diagram in Fig. [lO^a) shows in particular 
that an increase of the upstream flow U in presence of 
rain augments the erosion rate. Experiments with greater 
channel depth also give large erosion rates, but in the in- 
termediate regime, no clear trend is visible. 



9 



Weak Flow 



Strong Flow 



0.21 
0.15 



3 0.1 

O 



0.05 



0.2 



0.15 



1 0.1 



0.05 




• 










• 


• 


w 


(a) 




- • 

• • 

- • 

9 


• 

• 


> • 


o 

# 


O <~ 


) 

TO 

) 




Erosion Rate 




a 








• 




• 


(b) 




• • 


• 

• • 


• • 


O 




. ♦ 




•• 


• 


O : 


Channel Aspect Ratio 





than with Pi. Also increasing U appears to increase the 
aspect ratio, but it remains larger than most experiments 



u 



o 



0.02 0.04 0.06 

U 



FIG. 10: (Color online) (a) Parameter diagram comparing 
erosion rate r s of the different experiments processed in case 
of seepage erosion without rainfall (red hollow circle) and 
with rain (blue solid circle). The two parameters used are 
the strength of the upstream flow U = (H2 2 — H\ 2 )/L 2 and 
the size of the outlet divided by the length of experiment 
doutiet/L. The size of symbol is proportional to the erosion 
rate estimated from the volume measurement when the chan- 
nel length is 20 cm. The symbol scale size is ten times greater 
for measurements with rain, (b) Parameter diagram compar- 
ing channel aspect ratio A r (see text for definition) of the 
different experiments. The size of symbol is proportional to 
A r , and the scale is the same for experiments with P — and 
P = Pi- 



A channel aspect ratio can be defined as 
2max(C ?/ ) 



max(C ;E ) — mm(C x ) ' 



where, (C x , C y ) is a parametric representation in the con- 
tour d(x, y) = 1 cm. With this definition, A r = 1 for a 
semi-circle shaped channel and infinity for a line. Fig- 
ure [lO^b) displays A r for different values of experimental 
parameters when the channel length is equal to 25 cm. 
In order to compare A r across the various experiments, 
we evaluate it when channel length is equal to 25 cm and 
find that without local rain, the mean aspect ratio is 
< A r >= 2.04 and its standard deviation is OA r = 0.30, 
which means that for this channel length, the channel 
width is close to half of the length. Similarly, we find the 
values with rain to be < A r >= 1.54 and OA r = 0.22. 
Consequently channels in presence of homogeneous rain 
become wider, as expected. Although A r appears scat- 
tered in Fig. [lO^b), the cases with P = appear larger 



without rain. Overall, Fig. 10 appears to show that for 
P = Pi, the morphology of the channels created with 
a combination of rain and upstream flow are similar in 
shape to that with only rain. Altogether, these examples 
suggest that the overall channel evolution is not modi- 
fied by the presence of intermittent avalanches along the 
bank present in our experiments. 



VI. CONCLUSIONS 

Seepage erosion of a flat granular bed is experimentally 
studied at the laboratory scale to understand the influ- 
ence of the source of groundwater flow on the shape of 
channels. Significant differences are found between the 
case where the ground water comes primarily through 
a boundary from a far away source and the case where 
it is fed by uniform local rain. In the first case, when 
the ground water flow field should be governed by the 
Laplace equation, the channels formed are elongated with 
a roughly constant width and point towards the bound- 
ary through which water enters the bed. This shape 
could arise because the tip of the channel drains more wa- 
ter into the channel compared with the sides and there- 
fore most of the channel growth occurs at channel tip. By 
contrast, if uniform rain occurs on the granular bed, the 
seepage flow is expected to be governed now by the Pois- 
son equation and water is brought to the outlet almost 
uniformly from all directions in the half-plane. There- 
fore, the channel appears to grow uniformly roughly with 
a semi-circular shape. However, once the channel grows 
sufficiently long, local maxima of water flux entering the 
channel appear symmetrically on either side leading to 
a channel bifurcation. In addition, any perturbation of 
this growing front could in principle drain more water and 
therefore is unstable to a fingering instability in the con- 
text of the models using the Dupuit approximation |8j [9] . 
In case of shallow channels, we find the beginning of a 
splitting under rain which can be easily interpreted by 
this mechanism. This supports the notion that a chan- 
nel network can develop in a homogeneous bed whereby 
ground water flow splits as the channels grow leading the 
channels to split in turn. 

Further, we have shown the importance of the chan- 
nel depth in channelization, which is controlled by the 
parameter d out i e t' The upper part of the bed which is 
unsaturated becomes cohesive due to capillary bonds be- 
tween grains and collapse as a solid block as its founda- 
tions are sapped by the groundwater flow. But, while the 
spatio-temporal evolution of the erosion front appears to 
differ strongly, integrated quantities such as the erosion 
rate appear similar. 

Our observations have important implications for the 
interpretation of field data because numerous perturba- 
tions due to vegetation, pebbles, sediment inhomogeneity 
are present in nature that could influence channel dy- 
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namics. Perturbation of the erosion front due to random 
avalanching events are shown to not always lead to bifur- 
cations unless they are supported by underlying changes 
in ground water flow. This occurs because the seepage 
flow itself is not dependent on surface details. Therefore, 
our study shows that the overall shape of the channels 
is primarily determined by the nature of the groundwa- 
ter flow aided by the convergence of the flow near the 
channel tip. 

Our experiments also show significant effect of capillar- 
ity on the evolution of the channel bank which is relevant 
to beach rills and spring banks [I3j [14] . When the depth 
of the channel is smaller or close to the height of cap- 
illary rise, the bed is fully saturated and has a smooth 
appearance. For deeper channels, the upper part of the 
bed is unsaturated and cohesive due to capillary bonds 
between grains. In this case, the bank collapses as a solid 
block as its foundations are sapped by the groundwater 



flow. The capillary rise effects in the unsaturated region 
in our experiments may be also representative of similar 
effects in much larger sedimentary beds where clay can 
provide greater cohesion between grains and in channels 
observed in Navajo sandstones in Colorado [TO} [35]. 



Acknowledgments 

We thank Joshua Meyer who did preliminary mea- 
surements and Samuel King and Fouad Abdul Ameer 
for experimental assistance at the end of this project, 
and Daniel Abrams for stimulating discussions. This 
work was funded by the Department of Energy grants 
DE-FG0202ER15367(Clark) and DE-FG0299ER15004 
(MIT). 



[1] N. Izumi and G. Parker, J. Fluids Mech. 283, 341 (1995). 
[2] N. Izumi and G. Parker, J. Fluids Mech. 419, 239 (2000). 
[3] M. Reitz, D. Jerolmack, and J. Swenson, Geophysical 

Research Letters 37, L06401 (2010). 
[4] T. Dunne, Progress in Physical Geography 4, 211 (1980). 
[5] C. G. Higgins, Geology 10, 147 (1982). 
[6] J. Bear, Hydraulics of groundwater (McGraw-Hill New 

York, 1979). 

[7] S. A. Schumm, K. F. Boyd, C. G. Wolff, and W. Spitz, 

Geomorphology 12, 281 (1995). 
[8] D. M. Abrams, A. E. Lobkovsky, A. P. Petroff, K. M. 

Straub, B. McElroy, D. C. Mohrig, A. Kudrolli, and D. H. 

Rothman, Nature Geoscience 28, 193 (2009). 
[9] A. P. Petroff, O. Devauchelle, D. M. Abrams, A. E. 

Lobkovsky, A. Kudrolli, and D. H. Rothman, J. Fluids 

Mech. 673, 245 (2011). 
[10] M. P. Lamb, A. D. Howard, J. Johnson, K. X. Whipple, 

W. E. Dietrich, and J. T. Perron, J. Geophys. Res. Ill, 

E07002 (2006). 

[11] M. P. Lamb, W. E. Dietrich, S. M. Aciego, D. J. DePaolo, 
and M. Manga, Science 320, 1067 (2008). 

[12] P. D. Komar, Beach Processes and Sedimentation (Pren- 
tice Hall, 1998). 

[13] E. G. Otvos, Journal of Coastal Research 15, 1040 
(1999). 

[14] N. Schorghofer, B. Jensen, A. Kudrolli, and D. H. Roth- 
man, J. Fluids Mech. 503, 357 (2004). 

[15] A. D. Howard and C. F. McLane, Water Resources Re- 
search pp. 1659-1674 (1988). 

[16] A. D. Howard, Sapping Features of the Colorado Plateau: 
A Comparative Planetary Geology Field Guide pp. 71-83 
(1988). 

[17] B. Gomez and V. T. Mullen, Earth Surface Processes and 
Landforms 17, 465 (1992). 

[18] A. E. Lobkovsky, B. Jensen, A. Kudrolli, and D. H. Roth- 
man, J. Geophys. Res. 109, F04010 (2004). 

[19] A. E. Lobkovsky, B. E. Smith, A. Kudrolli, D. C. Mohrig, 



and D. H. Rothman, J. Geophys. Res. 112, F03S12 
(2007). 

[20] A. Pronprommin, Y. Takei, A. Mezgebu-Wubneh, and 
N. Izumi, Journal of hydro- environment Research 3, 232 
(2010). 

[21] A. Pronprommin and N. Izumi, J. Geophys. Res. 115, 
F02022 (2010). 

[22] A. Mezgebu-Wubneh, J. Nagahara, and N. Izumi, in 
River, Coastal and Estuarine Morphody namics 2011 (Ts- 
inghua University Press, Beijing, 2011), p. 368. 

[23] W.-J. Ni and H. Capart, J. Geophys. Res. Ill, F02014 
(2006). 

[24] A. E. Lobkovsky, A. Orpe, R. Molloy, A. E. Lobkovsky, 
A. Kudrolli, and D. H. Rothman, J. Fluids Mech. 605, 
47 (2008). 

[25] F. Charru and H. Mouilleron-Arnould, J. Fluids Mech. 

452, 303 (2002). 
[26] A. P. Petroff, O. Devauchelle, A. Kudrolli, and D. H. 

Rothman, Comptes Rendus Geoscience 344, 33 (2012). 
[27] S. Nowak, A. Samadani, and A. Kudrolli, Nature Physics 

1, 50 (2005). 

[28] L. Bocquet, E. Charlaix, and F. Restagno, Comptes 
Rendus de Physique de lAcademie des Sciences 3, 207 
(2002). 

[29] S. Herminghaus, Advances in Physics 54, 221 (2006). 

[30] G. A. Fox, M. L. Chu-Agor, and G. V. Wilson, Soil Sci- 
ence Society of America Journal 71, 1822 (2007). 

[31] S. Bonnet, Nature Geoscience 2, 766 (2009). 

[32] W. D. Ellison, Science 111, 245 (1950). 

[33] O. Devauchelle, A. P. Petroff, A. E. Lobkovsky, and D. H. 
Rothman, J. Fluids Mech. 667, 38 (2011). 

[34] M. Caponieri, S. Douady, S. Fauve, and C. Laroche, in 
Mobile Particulate Systems, edited by E. Guazzelli and 
L. Oger (Kluwer, Dordrecht, 1995), pp. 331-366. 

[35] M. P. Lamb, Ph.D. thesis, University of California, 
Berkeley (2008). 



